%%%%% Define constants
Tmax = 476;

%%%%% Add paths
addpath('../code/general_funcs');

%%%%% Load inputs
% Load consumer Adoptions
load('../data/data_mat/platformData/consumerAdoptionsData_CantonId_day.mat');
[Tall,K,~] = size(adoptionData.transactions);
initConsumers = adoptionData.pastConsumerTransactionsOrig'; % 1 x K
NewConsumers = reshape(sum(adoptionData.transactions,2), [Tall K]); % T x K
clear adoptionData;
consumersAtPeriodEnd = initConsumers + cumsum(NewConsumers,1);
consumersAtPeriodStart = consumersAtPeriodEnd - NewConsumers;

% Load provider movements
load('../data/data_mat/platformData/providerMovementsData_CantonId_day.mat');
NewProviders     = movementData.values(:,:,find(strcmp(movementData.VarNames,'NewProviders'))); % Tall x K
DeletedProviders = movementData.values(:,:,find(strcmp(movementData.VarNames,'DeletedProviders'))); % Tall x K
laggedProviders  = movementData.values(:,:,find(strcmp(movementData.VarNames,'LaggedProviders'))); % Tall x K
clear movementData;

initProviders = laggedProviders(1,:);
providersAtPeriodEnd = initProviders + cumsum(NewProviders-DeletedProviders,1);
assert(isequal(providersAtPeriodEnd-NewProviders+DeletedProviders, laggedProviders))

% Load covariates
load('../data/data_mat/covariates/covariates_CantonId_day.mat');
periodLabels = data.timeLabels;
geoLabels    = data.geoLabels;

myZnames = {'Area', 'Population', 'PopDensity', 'NumPersons19orOlder', 'NumHHsWithCar', 'NumTouristBeds', 'NumCommutersToLocation', ...
	'NumCarRentalAgencies', 'NumPublicTransitStops', 'MedianIncome', 'UnemploymentRate', 'PctageVoix_JOLY'};

[~,Zidxes] = ismember(myZnames, data.Znames);
Z = data.Z(:,Zidxes);
Znames = data.Znames(Zidxes);

%%%%%%% Compute density measures by area
idx = find(strcmp(Znames, 'Population'));             Population = Z(:,idx);
idx = find(strcmp(Znames, 'Area'));                   Area = Z(:,idx);
idx = find(strcmp(Znames, 'NumTouristBeds'));         TouristBedDensity      = Z(:,idx)./Area;
idx = find(strcmp(Znames, 'NumCommutersToLocation')); CommuterDensity        = Z(:,idx)./Area;
idx = find(strcmp(Znames, 'NumCarRentalAgencies'));   CarRentalAgencyDensity = Z(:,idx)./Area;
idx = find(strcmp(Znames, 'NumPublicTransitStops'));  TransitStopDensity     = Z(:,idx)./Area;
Z = [Z TouristBedDensity CommuterDensity CarRentalAgencyDensity TransitStopDensity];
Znames = [Znames {'TouristBedDensity', 'CommuterDensity', 'CarRentalAgencyDensity', 'TransitStopDensity'}];

% Process diploma level (aggregate)
belowBacDiplomaLevels = {'Fraction_DIPL_01','Fraction_DIPL_02','Fraction_DIPL_03','Fraction_DIPL_11','Fraction_DIPL_12','Fraction_DIPL_13','Fraction_DIPL_14'};
aboveBacDiplomaLevels = {'Fraction_DIPL_17','Fraction_DIPL_18'};
[~,below_dipl_idxes] = ismember(belowBacDiplomaLevels, data.Znames);
[~,above_dipl_idxes] = ismember(aboveBacDiplomaLevels, data.Znames);
FractionDiplBelowBac = sum(data.Z(:,below_dipl_idxes),2);
FractionDiplAboveBac = sum(data.Z(:,above_dipl_idxes),2);
Z      = [Z FractionDiplBelowBac FractionDiplAboveBac];
Znames = [Znames {'FractionDiplBelowBac', 'FractionDiplAboveBac'}];

% Base of participants at end of observation window
providersAtWindowEnd = providersAtPeriodEnd(Tmax,:)';
consumersAtWindowEnd = consumersAtPeriodEnd(Tmax,:)';
providerDensityAtWindowEnd = providersAtWindowEnd./Area;
consumerDensityAtWindowEnd = consumersAtWindowEnd./Area;
Z = [Z providersAtWindowEnd consumersAtWindowEnd providerDensityAtWindowEnd consumerDensityAtWindowEnd];
Znames = [Znames {'providersAtWindowEnd', 'consumersAtWindowEnd', 'providerDensityAtWindowEnd', 'consumerDensityAtWindowEnd'}];


%%%%% Make Table 3
T3 = summary_func(Z);
T3.Properties.RowNames = Znames;
T3 = T3({'Area', 'Population', 'PopDensity', 'NumPersons19orOlder', 'NumHHsWithCar', ...
	'providersAtWindowEnd', 'consumersAtWindowEnd', 'providerDensityAtWindowEnd', 'consumerDensityAtWindowEnd', ...
	'TouristBedDensity', 'CommuterDensity', 'CarRentalAgencyDensity', 'TransitStopDensity', ...
	'MedianIncome', 'UnemploymentRate', 'FractionDiplBelowBac', 'FractionDiplAboveBac', 'PctageVoix_JOLY'},[1 2 3 5 6]);
disp(T3);
writetable(T3, '../results/EDA/table3.csv', 'WriteRowNames', true);


%%%%% Make time series plots
aggConsumerBase  = sum(consumersAtPeriodEnd,2); % Tall x 1
aggProviderBase  = sum(providersAtPeriodEnd,2); % Tall x 1
aggNewConsumers     = sum(NewConsumers,2);           % Tall x 1
aggNewProviders     = sum(NewProviders,2);           % Tall x 1
aggDeletedProviders = sum(DeletedProviders,2);       % Tall x 1

%%% Fig 3
Mday = table(periodLabels, aggNewConsumers, aggNewProviders, aggDeletedProviders, aggConsumerBase, aggProviderBase);
plotTimeSeries(Mday(1:Tmax,[1 2:4]), '../results/EDA/fig3a.pdf', '', '', ...
    {'Consumer adoptions', 'Provider adoptions', 'Provider exits'}, {false false false}, {[0.35 0.35 0.35], [0,0,0], [0,0,0]}, {0.7,0.7,0.7}, 'northwest', [0, 230], '', {'-', '--', ':'}, 50, 20, 15);

plotTimeSeries(Mday(1:Tmax,[1 5:6]), '../results/EDA/fig3b.pdf', '', '', ...
    {'Base of consumers', 'Base of providers'}, {false false}, {[0.35 0.35 0.35], [0,0,0]}, {1,1}, 'northwest', [0, 20000], '', {'-', '--'}, 50, 20, 15);


%%%%% Fig W2
weekLabels = periodLabels(1:7:end);
aggConsumerBaseWeek = aggConsumerBase(7:7:end,:);
aggProviderBaseWeek = aggProviderBase(7:7:end,:);
NumWeeks = floor(size(aggConsumerBase,1)/7);
weekLabels = weekLabels(1:NumWeeks);
aggConsumerBaseWeek = aggConsumerBaseWeek(1:NumWeeks);
aggProviderBaseWeek = aggProviderBaseWeek(1:NumWeeks);
aggNewConsumersWeek     = zeros(NumWeeks,1);
aggNewProvidersWeek     = zeros(NumWeeks,1);
aggDeletedProvidersWeek = zeros(NumWeeks,1);
for ww = 1:NumWeeks
	ww2_dds = (ww-1)*7+1:(ww-1)*7+7; % Days corresponding to week ww
	aggNewConsumersWeek(ww)     = sum(aggNewConsumers(ww2_dds));
	aggNewProvidersWeek(ww)     = sum(aggNewProviders(ww2_dds));
	aggDeletedProvidersWeek(ww) = sum(aggDeletedProviders(ww2_dds));
end

isWindowLimit = zeros(NumWeeks,1); isWindowLimit(68) = 1;

% Add some checks
assert(isequal(aggConsumerBaseWeek(2:end), aggConsumerBaseWeek(1) + cumsum(aggNewConsumersWeek(2:end))));
aggProvidersDeltaWeek = aggNewProvidersWeek - aggDeletedProvidersWeek;
assert(isequal(aggProviderBaseWeek(2:end), aggProviderBaseWeek(1) + cumsum(aggProvidersDeltaWeek(2:end))));

Mweek = table(weekLabels, aggNewConsumersWeek, aggNewProvidersWeek, aggDeletedProvidersWeek, aggConsumerBaseWeek, aggProviderBaseWeek, isWindowLimit);

%%%% DEFINE TICKS FOR WEEK-LEVEL TIME SERIES OVER ENTIRE DATA WINDOW
xticks_values = [1 27 53 79 105 131 157 183 209 235]';
xticks_values_disp ={'2013-01-01', '2013-07-01', '2014-01-01', '2014-07-01', '2015-01-01', '2015-07-01', ...
	'2016-01-01', '2016-07-01', '2017-01-01', '2017-07-01'}';

tt_windowLimit = [68];

verticalLineColor = [0.6,0.6,0.6];
verticalLineWidth = 3;


%%% Fig W2a
values = table2array(Mweek(:,[2:4]));
NumVars = size(values,2);
varNames4legend = {'Consumer adoptions', 'Provider adoptions', 'Provider exits'};
legendLocation = 'northwest';
ylimits = [0 3000];
lineWidths = {0.7, 0.7, 0.7};
lineStyles = {'-', '--', ':'};
lineColors = {[0.35 0.35 0.35], [0,0,0], [0,0,0]};
axisTitleFontSize = 20;
xtickLabelFontSize = 15;

% Make figure + plot time series and vertical indicators
h = figure;

set(gcf, 'Units', 'Inches', 'Position', [0, 0, 10, 7], 'PaperUnits', 'Inches', 'PaperSize', [10, 7]);
for vv = 1:NumVars
	myvalues = values(:,vv);
 	plot(myvalues,'LineWidth', lineWidths{vv}, 'LineStyle', lineStyles{vv}, 'Color', lineColors{vv});;
	ylim(ylimits);
	hold on;
end
%%% Add vertical line for end of period of interest
ylimits = ylim;
ymin = ylimits(1);
ymax = ylimits(2);
plot([tt_windowLimit tt_windowLimit],[ymin ymax], 'LineWidth', verticalLineWidth, 'Color', verticalLineColor, 'LineStyle', ':');
hold off;

xlim([0,NumWeeks]);

% Make legend	
[legh,objh] = legend(varNames4legend, 'Location', legendLocation, 'FontSize',16);
	
% Print xticks
xticks(xticks_values); xticklabels(xticks_values_disp); xtickangle(45);
a = get(gca,'XTickLabel');
set(gca,'XTickLabel',a,'fontsize',xtickLabelFontSize)

% Reformat yticks to avoid scientific notation
set(gca,'YTickLabel',formatNumbersWith1000commas(get(gca,'YTick'))');	


%% Add annotation
% Arrow
annotation('doublearrow', [0.13055,0.3333], [0.6,0.6]);
% Textbox
mytextbox = annotation('textbox', 'String', 'Time window analyzed', 'Position', [0.165 0.5 0.15 0.2], 'FontSize', 16, 'EdgeColor', 'none', 'FitBoxToText','off');

% Save to output file
print(h, '-dpdf', '../results/EDA/figW2a.pdf');
close(h);


%%% Fig W2b
values = table2array(Mweek(:,[5:6]));
NumVars = size(values,2);
varNames4legend = {'Base of consumers', 'Base of providers'};
legendLocation = 'northwest';
ylimits = [0 250000];
lineWidths = {1, 1};
lineStyles = {'-', '--'};
lineColors = {[0.35 0.35 0.35], [0,0,0]};
axisTitleFontSize = 20;
xtickLabelFontSize = 15;

% Make figure + plot time series and vertical indicators
h = figure;

set(gcf, 'Units', 'Inches', 'Position', [0, 0, 10, 7], 'PaperUnits', 'Inches', 'PaperSize', [10, 7]);
for vv = 1:NumVars
	myvalues = values(:,vv);
 	plot(myvalues,'LineWidth', lineWidths{vv}, 'LineStyle', lineStyles{vv}, 'Color', lineColors{vv});;
	ylim(ylimits);
	hold on;
end
%%% Add vertical line for end of period of interest
ylimits = ylim;
ymin = ylimits(1);
ymax = ylimits(2);
plot([tt_windowLimit tt_windowLimit],[ymin ymax], 'LineWidth', verticalLineWidth, 'Color', verticalLineColor, 'LineStyle', ':');
hold off;

xlim([0,NumWeeks]);

% Make legend	
[legh,objh] = legend(varNames4legend, 'Location', legendLocation, 'FontSize',16);
	
% Print xticks
xticks(xticks_values); xticklabels(xticks_values_disp); xtickangle(45);
a = get(gca,'XTickLabel');
set(gca,'XTickLabel',a,'fontsize',xtickLabelFontSize)

% Reformat yticks to avoid scientific notation
set(gca,'YTickLabel',formatNumbersWith1000commas(get(gca,'YTick'))');	


%% Add annotation
% Arrow
annotation('doublearrow', [0.13055,0.3333], [0.6,0.6]);
% Textbox
mytextbox = annotation('textbox', 'String', 'Time window analyzed', 'Position', [0.165 0.5 0.15 0.2], 'FontSize', 16, 'EdgeColor', 'none', 'FitBoxToText','off');

% Save to output file
print(h, '-dpdf', '../results/EDA/figW2b.pdf');
close(h);



%%%%% Prepare everything to make maps for Fig 4
mydayLabels = {'2013-01-01', '2013-09-01', '2014-04-01'};
[~,mydayIdxes] = ismember(mydayLabels, periodLabels);

installedBases = zeros(K,length(mydayIdxes),2);
installedBases(:,:,1) = consumersAtPeriodEnd(mydayIdxes,:)'; % K x 3
installedBases(:,:,2) = providersAtPeriodEnd(mydayIdxes,:)'; % K x 3

for ll = 1:2
	if ll == 1; userType = 'consumer'; varLabel = 'CumulativeConsumersPer100000People'; end
	if ll == 2; userType = 'provider'; varLabel = 'CumulativeProvidersPer100000People'; end
	for tidx = 1:length(mydayIdxes)
		tt_label = mydayLabels{tidx};
		outputFile = sprintf('../results/EDA/adoptionMaps/InstalledBasePerCapita_%s_%s.csv', userType, tt_label);
		myTable = table(geoLabels, 1e5*installedBases(:,tidx,ll)./Population);
		myTable.Properties.VariableNames = {'CantonId', varLabel};
		writetable(myTable, outputFile);
	end
end
